t 2 g $$ {t}_{2g} $$ d $$ d $$ orbital ordering patterns in KBF 3 $$ {}_3 $$ (B = Sc, Ti, Fe, Co) perovskites.

The orbital ordering (OO) resulting from the partial occupancy of the t 2 g $$ {t}_{2g} $$ d $$ d $$ subshell of the transition metals in KBF 3 $$ {}_3 $$ (B = Sc, Ti, Ffe, Co) perovskites, and the many possible patterns arising from the coupling between the B sites, have been investigated at the quantum mechanical level ( all electron $$ \mathrm{all}\ \mathrm{electron} $$ Gaussian type basis set, B3LYP hybrid functional) in a 40 atoms supercell. The numerous patterns are distributed into 162 classes of equivalent configurations. For each fluoroperovskite, one representative per class has been calculated. The four compounds behave similarly: an identical dependence of the energy and volume (or cell parameters) on the OO pattern; the spanned energy interval is small (1 to 2 mE h $$ {}_h $$ per formula unit), suggesting that most of the configurations are occupied at room and even at low temperature. A linear model, taking into account the relative orbital order in contiguous sites, reproduces the energy order in the full set for each compound, suggesting that it could be used for studying OO in larger supercells.

The energies and charge and spin distributions in the low-lying levels of singlet and triplet N2V defects in diamond from direct variational calculations of the excited states.

This paper reports the energies and charge and spin distributions of the low-lying excited states in singlet and triplet N2V defects in diamond from direct Δ-SCF calculations based on Gaussian orbitals within the B3LYP, PBE0, and HSE06 functionals. They assign the observed absorption at 2.463 eV, first reported by Davies et al. [Proc. R. Soc. London 351, 245 (1976)], to the excitation of a N(sp3) lone-pair electron in the singlet and triplet states, respectively, with estimates of ∼1.1 eV for that of the unpaired electrons, C(sp3). In both cases, the excited states are predicted to be highly local and strongly excitonic with 81% of the C(sp3) and 87% of the N(sp3) excited charges localized at the three C atoms nearest neighbor (nn) to the excitation sites. Also reported are the higher excited gap states of both the N lone pair and C unpaired electron. Calculated excitation energies of the bonding sp3 hybrids of the C atoms nn to the four inner atoms are close to that of the bulk, which indicates that the N2V defect is largely a local defect. The present results are in broad agreement with those reported by Udvarhelyi et al. [Phys. Rev. B 96, 155211 (2017)] from plane wave HSE06 calculations, notably for the N lone pair excitation energy, for which both predict an energy of ∼2.7 eV but with a difference of ∼0.5 eV for the excitation of the unpaired electron.

Catalytic CO2 Reduction with Heptacoordinated Polypyridine Complexes: Switching the Selectivity via Metal Replacement.

The discovery of molecular catalysts for the CO2 reduction reaction (CO2 RR) in the presence of water, which are both effective and selective towards the generation of carbon-based products, is a critical task. Herein we report the catalytic activity towards the CO2 RR in acetonitrile/water mixtures by a cobalt complex and its iron analog both featuring the same redox-active ligand and an unusual seven-coordination environment. Bulk electrolysis experiments show that the cobalt complex mainly yields formate (52 % selectivity at an applied potential of -2.0 V vs Fc+ /Fc and 1 % H2 O) or H2 (up to 86 % selectivity at higher applied bias and water content), while the iron complex always delivers CO as the major product (selectivity >74 %). The different catalytic behavior is further confirmed under photochemical conditions with the [Ru(bpy)3 ]2+ sensitizer (bpy=2,2'-bipyridine) and N,N-diisopropylethylamine as electron donor, where the cobalt complex leads to preferential H2 formation (up to 89 % selectivity), while the iron analog quantitatively generates CO (up to 88 % selectivity). This is ascribed to a preference towards a metal-hydride vs. a metal-carboxyl pathway for the cobalt and the iron complex, respectively, and highlights how metal replacement may effectively impact on the reactivity of transition metal complexes towards solar fuel formation.

The role of the A monovalent cation in the AVF3 perovskite series. A quantum mechanical investigation.

The relative stability of various phases of five AVF3 compounds (A = Li, Na, K, Rb and Cs) is investigated starting from the cubic (C) Pm3̄m (221) prototype structure, with five atoms (one formula unit) in the primitive cell. To the authors' knowledge, only three of these compounds have been investigated experimentally (Na, K and Rb), and they are reported as being cubic. The picture emerging from the present simulation is quite different: CsVF3 and RbVF3 are dynamically stable in the cubic structure, KVF3 is tetragonal, with space group (SG) I4/mcm (no. 140) and 10 atoms in the unit cell; note, however, that an orthorhombic Pnma (62) phase (four formula units) exists, which is not a subgroup of I4/mcm (140), and is very close in energy to the tetragonal phase. Further symmetry lowering is observed in the Na and Li compounds that are orthorhombic. The energy gain and volume reduction with respect to the cubic aristotype increase along the series K, Na and Li, and is very large for the last two compounds. Both FM and AFM solutions have been explored, and they show a very similar path along the SG modifications. The present scheme for determining the lowest energy SG is general, and can be applied to any perovskite. The B3LYP full range hybrid functional and the Hartree-Fock (HF) Hamiltonian, an all-electron Gaussian type basis set and the CRYSTAL code have been used.

First-principles study of the structural and electronic properties of tetragonal ZrOX (X = S, Se, and Te) monolayers and their vdW heterostructures for applications in optoelectronics and photocatalysis.

Using first-principles calculations, we have studied the structural and electronic properties of ZrOX (X = S, Se, and Te) monolayers and their van der Waals heterostructures in the tetragonal structure. Our results show that these monolayers are dynamically stable and are semiconductors with electronic bandgaps ranging from 1.98 to 3.16 eV as obtained with the GW approximation. By computing their band edges, we show that ZrOS and ZrOSe are of interest for water splitting applications. In addition, the van der Waals heterostructures formed by these monolayers show a type I band alignment for ZrOTe/ZrOSe and a type II alignment for the other two heterostructures, making them potential candidates for certain optoelectronic applications involving electron/hole separation.

The d Orbital Multi Pattern Occupancy in a Partially Filled d Shell: The KFeF3 Perovskite as a Test Case.

The occupancy of the d shell in KFeF3 is t2g4eg2, with five α and one ß electrons. The Jahn-Teller lift of degeneracy in the t2g sub-shell produces a tetragonal relaxation of the unit cell (4.09 vs. 4.22 Å, B3LYP result) not observed experimentally. In order to understand the origin of this apparent contradiction, we explored, with a 2 × 2 × 2 supercell (40 atoms per cell), all possible local structures in which contiguous Fe atoms have a different occupancy of the t2g orbitals with the minority spin electron. A total of 6561 configurations (with occupancies from (8,0,0) to (3,2,2) of the 3 t2g orbitals of the 8 Fe atoms) have been explored, with energies in many cases lower (by up to 1550 µEh per 2 Fe atoms) than the one of the fully ordered case, both for the ferromagnetic and the anti-ferromagnetic solutions. The results confirm that the orientation of the ß d electron of Fe influences the electrostatics (more efficient relative orientation of the Fe quadrupoles of the d shell) of the system, but not the magnetic interactions. Three hybrid functionals, B3LYP, PBE0, and HSE06, provide very similar results.

How deeply are core electrons perturbed when valence electrons are spin polarized? The case study of transition metal compounds.

The ferromagnetic and antiferromagnetic wave functions of the KMnF3 perovskite have been evaluated quantum-mechanically by using an all electron approach and, for comparison, pseudopotentials on the transition metal and the fluorine ions. It is shown that the different number of α and ß electrons in the d shell of Mn perturbs the inner shells, with shifts between the α and ß eigenvalues that can be as large as 6 eV for the 3s level, and is far from negligible also for the 2s and 2p states. The valence electrons of F are polarized by the majority spin electrons of Mn, and in turn, spin polarize their 1s electrons. When a pseudopotential is used, such a spin polarization of the core functions of Mn and F can obviously not take place. The importance of such a spin polarization can be appreciated by comparing (i) the spin density at the Mn and F nuclear position, and then the Fermi contact constant, a crucial quantity for the hyperfine coupling, and (ii) the ferromagnetic-antiferromagnetic energy difference, when obtained with an all electron or a pseudopotential scheme, and exploring how the latter varies with pressure. This difference is as large as 50% of the all electron datum, and is mainly due to the rigid treatment of the F ion core. The effect of five different functionals on the core spin polarization is documented.

The role of spin density for understanding the superexchange mechanism in transition metal ionic compounds. The case of KMF3 (M = Mn, Fe, Co, Ni, Cu) perovskites.

In many recent papers devoted to first row transition metal fluorides and oxides, not much attention is devoted to the spin density, a crucial quantity for the determination of the superexchange mechanism, and then for the ferro-antiferromagnetic energy difference. Usually, only the eigenvalues of the system are represented, in the form of band structures or, more frequently, of density of states (DOS). When discussing the orbital ordering and the Jahn-Teller effect, simple schemes with cubes and lobes are used to illustrate the shape of the d occupancy. But the eigenvectors, and the resulting spin density function, as obtained from the calculations, are rarely shown. When represented, only a fuzzy shape that recalls the d orbital shape can be observed. On the basis of these considerations, spin density maps for 5 compounds of the KMF3 (M from Mn to Cu) family have been produced, which clearly illustrate which d orbital is singly or doubly occupied. At variance with respect to the near totality of the papers devoted to these systems, we use an all electron scheme, a Gaussian type basis set, and the Hartree-Fock Hamiltonian or the B3LYP hybrid functional (the resulting maps turn out to be very similar, in the scale used for our figures). The spin density in the five cases can easily be interpreted in terms of the shape of the d orbitals as appearing in textbooks.

Quantum mechanical simulation of various phases of KVF3perovskite.

The relative stability ΔEof the cubicPm3¯m(C), of the two tetragonalP4mbm(T1) andI4mcm(T2), and of the orthorhombicPbnm(O) phases of KVF3has been computed both for the ferromagnetic (FM) and antiferromagnetic (AFM) solutions, by using the B3LYP full range hybrid functional and the Hartree-Fock (HF) Hamiltonian, an all-electron Gaussian type basis set and the CRYSTAL code. The stabilization of the T2 phase with respect to the C one (152µHa for B3LYP, 180µHa for HF, per 2 formula units) is due to the rotation of the VF6octahedra with respect to thecaxis, by 4.1-4.6 degrees. The O phase is slightly less stable than the T2 phase (by 6 and 20µHa for B3LYP and HF); it is, however, a stable structure as the dynamical analysis confirms. The mechanism of the stabilization of the AFM solution with respect to the FM one is discussed through the spin density maps, and is related to the key role of the exact exchange term (20% in B3LYP, 100% in HF). The G-AFM phase (the first six neighbors of the reference V ion with spin reversed) is more stable than the FM one by about 500 (HF) and 1800 (B3LYP)µHa per two formula units. A volume reduction is observed in the C to T passage, and in the FM to AFM one, both being of the order of 0.3-0.5A˚3at the B3LYP level. Atomic charges, magnetic moments and bond populations, evaluated according to a Mulliken partition of the charge a spin density functions, complete the analysis. The IR and Raman spectra of the FM and AFM C, T2 and O cells are discussed; the only noticeable difference between the various space groups appears in the modes with wavenumbers lower than 100 cm-1.

Assessing the Accuracy of Machine Learning Thermodynamic Perturbation Theory: Density Functional Theory and Beyond.

Machine learning thermodynamic perturbation theory (MLPT) is a promising approach to compute finite temperature properties when the goal is to compare several different levels of ab initio theory and/or to apply highly expensive computational methods. Indeed, starting from a production molecular dynamics trajectory, this method can estimate properties at one or more target levels of theory from only a small number of additional fixed-geometry calculations, which are used to train a machine learning model. However, as MLPT is based on thermodynamic perturbation theory (TPT), inaccuracies might arise when the starting point trajectory samples a configurational space which has a small overlap with that of the target approximations of interest. By considering case studies of molecules adsorbed in zeolites and several different density functional theory approximations, in this work we assess the accuracy of MLPT for ensemble total energies and enthalpies of adsorption. It is shown that problematic cases can be detected even without knowing reference results and that even in these situations it is possible to recover target level results within chemical accuracy by applying a machine-learning-based Monte Carlo (MLMC) resampling. Finally, on the basis of the ideas developed in this work, we assess and confirm the accuracy of recently published MLPT-based enthalpies of adsorption at the random phase approximation level, whose high computational cost would completely hinder a direct molecular dynamics simulation.

Strategies for the optimization of the structure of crystalline compounds.

When different proposals exist (or can reasonably be formulated) for the size of the unit cell (in terms of number of atoms) and space group of crystalline compounds, a strategy for exploring with simulation methods the various cases and for investigating their relative stability must be defined. The optimization schemes of periodic quantum mechanical codes work in fact at fixed space group and number of atoms per unit cell, so that only the fractional coordinates of the atoms and the lattice parameters are optimized. A strategy is here presented, based on four standard tools, used synergistically and in sequence: (1) the optimization of inner coordinates and unit cell parameters; (2) the calculation of the vibrational frequencies not only at Γ , but also at a set of k → points (in the example presented here they are eight, generated by a shrinking factor 2), looking for possible negative wavenumbers. The latter correspond to maxima, rather than minima, along the coordinate described by the corresponding normal mode; (3) the exploration of the total energy along the mode with negative wavenumber, looking for the minimum of the curve; (4) the identification of the new space group corresponding to the reduced symmetry resulting from the previous step. The strategy is illustrated with reference to the KMnF3 perovskite compound, for which many space groups are proposed in the literature, ranging from cubic Pm 3 ¯ m to tetragonal P 4 m bm or I 4 m cm and orthorhombic (Pnma and Cmcm) down to monoclinic (P21 /m). The corresponding primitive cells contain 5, 10, and 20 atoms in the various cases, and the point symmetry reduces from 48 to 4 operators. In nature, the KMnF3 phase transitions also include the magnetic phases. For simplicity, here we limit the analysis to the ones that take place between ferromagnetic phases, as they are sufficiently rich for illustrating the proposed strategy. As the total energy differences involved can be as small as, say, 10-50 µHartree, a high numerical accuracy at each one of the steps mentioned above is required. The present calculations, performed with the CRYSTAL code, by using an all electron basis set and the Hartree-Fock and B3LYP functionals, document such an accuracy. The energy difference between the tetragonal I 4 m cm and cubic Pm 3 ¯ m phases is 225 µHartree, with a volume reduction of 0.58 Å3 ; the differences between the orthorhombic and tetragonal phases are an order of magnitude smaller, being 23 µHartree and 0.06 Å3 for total energy and cell volume, respectively.

The ferromagnetic and anti-ferromagnetic phases (cubic, tetragonal, orthorhombic) of KMnF3. A quantum mechanical investigation.

Many space groups are proposed in the literature for the KMnF3 perovskite (see, for example, Knight et al., J. Alloys Compd., 2020, 842, 155935), ranging from cubic (C) (Pm3̄m) to tetragonal (T) ( or I4/m) down to orthorhombic (O) (Pbnm). The relative stability ΔE of these phases, both ferromagnetic (FM) and antiferromagnetic (AFM), has been investigated quantum mechanically by using both the B3LYP hybrid functional and the Hartree-Fock Hamiltonian, an all-electron Gaussian type basis set and the CRYSTAL code. The O phase is slightly more stable than the T phase which in turn is more stable than the C phase, in agreement with experimental evidence. The C to T to O transition is accompanied by a volume reduction. The mechanism of stabilization of the AFM solution with respect to the FM one is discussed. Spin density maps and profiles, Mulliken charges, magnetic moments and bond population data are used for supporting the proposed mechanism. The IR and Raman spectra of the FM and AFM C, T and O cells are discussed; the only noticeable difference between the C, T and O spectra appears at wavenumbers lower than 150 cm-1. The effect of pressure is also explored in the 0-20 GPa interval. The stability order (O > T > C) at 0 GPa persists also at high pressure, and the differences between the phases increase.

The superexchange mechanism in crystalline compounds. The case of KMF3(M = Mn, Fe, Co, Ni) perovskites.

The ferromagnetic and antiferromagnetic wavefunctions of four KMF3(M = Mn, Fe, Co and Ni) perovskites have been obtained quantum-mechanically with the CRYSTAL code, by using the Hartree-Fock (HF) Hamiltonian and three flavours of DFT (PBE, B3LYP and PBE0) and anall-electronGaussian type basis set. In the Fe and Co cases, with d6and d7occupation, the Jahn-Teller distortion of the cubic cell is as large as 0.12 Å. Various features of the superexchange interaction energies (SIE), namely additivity, dependence on the M-M distance, on theMFM̂angle, and on the adopted functional, are explored. The contribution to SIE by the Coulomb, exchange and kinetic energy terms is analyzed. It is shown that, when using density functionals, SIE clearly correlates with the amount of exact (Hartree-Fock) exchange in the functional. The effect of SIE on the equilibrium geometry and volume of the unit cell is discussed, and it is shown that the key quantity is the spin polarization of the (closed shell) F ions along the M-F-M path. The effect of thismagneticpressureis evaluated quantitatively for the first time. The superexchange coupling constantJ, evaluated at the HF level and through the Ising model, underestimates the experimental values by about 60%-70%. The more sophisticated Yamaguchi model (that takes into account the contamination of the FM and AFM spin states) does not reduce the discrepancy. The B3LYP hybrid functional overestimates the experiments. These last are bracketed by HF and PBE0. For PBE, the overestimation is huge. Finally, Mulliken population data, charge and spin density maps and density of states are used to illustrate the electronic structure.

Oxygen and vacancy defects in silicon. A quantum mechanical characterization through the IR and Raman spectra.

The Infrared (IR) and Raman spectra of various defects in silicon, containing both oxygen atoms (in the interstitial position, Oi) and a vacancy, are computed at the quantum mechanical level by using a periodic supercell approach based on a hybrid functional (B3LYP), an all-electron Gaussian-type basis set, and the Crystal code. The first of these defects is VO: the oxygen atom, twofold coordinated, saturates the unpaired electrons of two of the four carbon atoms on first neighbors of the vacancy. The two remaining unpaired electrons on the first neighbors of the vacancy can combine to give a triplet (Sz = 1) or a singlet (Sz = 0) state; both states are investigated for the neutral form of the defect, together with the doublet solution, the ground state of the negatively charged defect. Defects containing two, three, and four oxygen atoms, in conjunction with the vacancy V, are also investigated as reported in many experimental papers: VO2 and VOOi (two oxygen atoms inside the vacancy, or one in the vacancy and one in interstitial position between two Si atoms) and VO2Oi and VO22Oi (containing three and four oxygen atoms). This study integrates and complements a recent investigation referring to Oi defects [Gentile et al., J. Chem. Phys. 152, 054502 (2020)]. A general good agreement is observed between the simulated IR spectra and experimental observations referring to VOx (x = 1-4) defects.

Hybrid localized graph kernel for machine learning energy-related properties of molecules and solids.

Nowadays, the coupling of electronic structure and machine learning techniques serves as a powerful tool to predict chemical and physical properties of a broad range of systems. With the aim of improving the accuracy of predictions, a large number of representations for molecules and solids for machine learning applications has been developed. In this work we propose a novel descriptor based on the notion of molecular graph. While graphs are largely employed in classification problems in cheminformatics or bioinformatics, they are not often used in regression problem, especially of energy-related properties. Our method is based on a local decomposition of atomic environments and on the hybridization of two kernel functions: a graph kernel contribution that describes the chemical pattern and a Coulomb label contribution that encodes finer details of the local geometry. The accuracy of this new kernel method in energy predictions of molecular and condensed phase systems is demonstrated by considering the popular QM7 and BA10 datasets. These examples show that the hybrid localized graph kernel outperforms traditional approaches such as, for example, the smooth overlap of atomic positions and the Coulomb matrices.

Interstitial carbon defects in silicon. A quantum mechanical characterization through the infrared and Raman spectra.

The Infrared (IR) and Raman spectra of various interstitial carbon defects in silicon are computed at the quantum mechanical level by using an all electron Gaussian type basis set, the hybrid B3LYP functional and the supercell approach, as implemented in the CRYSTAL code (Dovesi et al. J. Chem. Phys. 2020, 152, 204111). The list includes two 〈100〉 split interstitial IXY defects, namely ICC and ICSi , a couple of related defects that we indicate as IX IY , the so called C i C s 0 in its A and B form, as well as SiCi Si and Cs Ci Cs , in which the interstitial carbon atom is twofold coordinated. The second undergoes a large relaxation, and the final configuration is close to ICC Cs . Geometries, relative stabilities, electronic, and vibrational properties are analysed. All these defects show characteristic features in their IR spectrum (above 730 cm- 1 ), whereas the Raman spectrum is dominated, in most of the cases, by the pristine silicon peak at 530 cm-1 , that hides the defect peaks.

The CRYSTAL code, 1976-2020 and beyond, a long story.

CRYSTAL is a periodic ab initio code that uses a Gaussian-type basis set to express crystalline orbitals (i.e., Bloch functions). The use of atom-centered basis functions allows treating 3D (crystals), 2D (slabs), 1D (polymers), and 0D (molecules) systems on the same grounds. In turn, all-electron calculations are inherently permitted along with pseudopotential strategies. A variety of density functionals are implemented, including global and range-separated hybrids of various natures and, as an extreme case, Hartree-Fock (HF). The cost for HF or hybrids is only about 3-5 times higher than when using the local density approximation or the generalized gradient approximation. Symmetry is fully exploited at all steps of the calculation. Many tools are available to modify the structure as given in input and simplify the construction of complicated objects, such as slabs, nanotubes, molecules, and clusters. Many tensorial properties can be evaluated by using a single input keyword: elastic, piezoelectric, photoelastic, dielectric, first and second hyperpolarizabilities, etc. The calculation of infrared and Raman spectra is available, and the intensities are computed analytically. Automated tools are available for the generation of the relevant configurations of solid solutions and/or disordered systems. Three versions of the code exist: serial, parallel, and massive-parallel. In the second one, the most relevant matrices are duplicated on each core, whereas in the third one, the Fock matrix is distributed for diagonalization. All the relevant vectors are dynamically allocated and deallocated after use, making the code very agile. CRYSTAL can be used efficiently on high performance computing machines up to thousands of cores.

The spectroscopic characterization of interstitial oxygen in bulk silicon: A quantum mechanical simulation.

The vibrational Infrared and Raman spectra of six interstitial oxygen defects in silicon containing a Si-O-Si bridge between adjacent Si atoms are obtained from all-electron B3LYP calculations within a supercell scheme, as embodied in the CRYSTAL code. Two series of defects have been considered, starting from the single interstitial defect, O1. The first consists of four defects, O1,n, in which two O1 defects are separated by (n - 1) Si atoms, up to n = 4. The second consists of four defects, On, in which nO1 defects surround a single Si atom, with n = 1-4, where O4 has the same local nearest neighbor structure as α-quartz. For both series of defects, the equilibrium geometries, charge distributions, and band structures are reported and analyzed. The addition of 1-4 oxygen atoms to the perfect lattice generates 3-12 new vibrational modes, which, as a result of the lighter atomic mass of O with respect to Si, are expected to occur at wavenumbers higher than 521 cm-1, the highest frequency of pristine silicon, thereby generating a unique new Raman spectrum. However, only a small subset of these new modes is found in the spectrum. They appear at 1153 cm-1 (O1), at 1049 cm-1 and 1100 cm-1 (O1,2), at 1108 cm-1 (O1,3), at 1130 cm-1 and 1138 cm-1 (O1,4), and 773 cm-1, 1057 cm-1, and 1086 cm-1 (O4), and can be considered "fingerprints" of the respective defects, as they are sufficiently well separated from each other. Graphical animations indicate the nature and intensity of each of the observed modes which are not overtones or combinations.

Nitrogen substitutional defects in silicon. A quantum mechanical investigation of the structural, electronic and vibrational properties.

The vibrational infrared (IR) and Raman spectra of seven substitutional defects in bulk silicon are computed, by using the quantum mechanical CRYSTAL code, the supercell scheme, an all electron Gaussian type basis set and the B3LYP functional. The relative stability of various spin states has been evaluated, the geometry optimized, the electronic structure analyzed. The IR and Raman intensities have been evaluated analitically. In all cases the IR spectrum is dominated by a single N peak (or by two or three peaks with very close wavenumbers), whose intensity is at least 20 times larger than the one of any other peak. These peaks fall in the 645-712 cm-1 interval, and a shift of few cm-1 is observed from case to case. The Raman spectrum of all defects is dominated by an extremely intense peak at about 530 cm-1, resulting from the (weak) perturbation of the peak of pristine silicon.

On the Models for the Investigation of Charged Defects in Solids: The Case of the VN- Defect in Diamond.

Local charged defects in periodic systems are usually investigated by adopting the supercell charge compensated (CC) model, which consists of two main ingredients: (i) the periodic supercell, hopefully large enough to reduce to negligible values the interaction among defects belonging to different cells; (ii) a background of uniform compensating charge that restores the neutrality of the supercell and then avoids the "Coulomb catastrophe". Here, an alternative approach is proposed and compared to CC, the double defect (DD) model, in which another point defect is introduced in the supercell that provides (or accept) the electron to be transferred (subtracted) to the defect of interest. The DD model requires obviously a (much) larger supercell than CC, and the effect of the relative position of the two defects must be explored. A third possible option, the cluster approach, is not discussed here. The two models have been compared with reference to the VN- defect; for DD, the positive compensating charge is provided by a P atom. Three cubic supercells of increasing size (containing 216, 512, and 1000 atoms) and up to eight relative VN--P+ defect-defect positions have been considered. The comparison extends to the equilibrium geometry around the defect, band structure, charge and spin distribution, IR and Raman vibrational spectra, and electron paramagnetic resonance constants. It turns out that the CC and DD models provide very similar results for all of these properties, in particular when the P+ compensating defect is not too close to VN-.